Classical field techniques for condensates in one-dimensional rings at finite temperatures 
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For a condensate in a one-dimensional ring geometry, we compare the thermodynamic properties of three 
conceptually different classical field techniques: stochastic dynamics, microcanonical molecular dynamics, and 
the "classical field method." Starting from non-equilibrium initial conditions, all three methods approach steady 
states whose distribution and correlation functions are in excellent agreement with an exact evaluation of the par- 
tition function in the high-temperature limit. Our study helps to establish these various classical field techniques 
as powerful non-perturbative tools for systems at finite temperatures. 
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I. INTRODUCTION 

Future experiments on ultracold gases will require theoret- 
ical input from numerical simulations which are reliable in 
the case of strong interactions beyond the realm of perturba- 
tive expansions. Powerful non-perturbative techniques have 
been developed ranging from the partial summation of Feyn- 
man diagrams in the context of condensed matter and nu- 
clear theory (e.g. the many-body T-matrix approach |l|]) to 
more recent renormalization group ideas inspired from quan- 
tum information theory (e.g. the density-matrix renormaliza- 
tion group (DMRG) 131). These methods, however, are of- 
ten intractable in systems at finite temperatures or beyond one 
spatial dimension. 

Many of the complications at finite temperatures arise from 
a need to account for higher-order classical fluctuations, for 
instance, near a second-order phase transition. Here, a rather 
broad category of methods, which we will refer to as "clas- 
sical field" (CF) methods, have proved to be a useful alter- 
native. Initially conceived of as a convenient way to sample 
the thermodynamic partition function, similar to Monte-Carlo 
techniques [3], a variant of these CF methods, known as the 
"classical field method" flrjll . has recently been developed and 
employed to study both the thermodynamics and real-time dy- 
namics of cold atomic systems. Some recent examples are 
studies of vortex-antivortex binding in the Kosterlitz-Thouless 
transition J3l or the shift of the Bose-Einstein condensation 
(BEC) transition temperature in an interacting Bose gas [5]. 
Despite its relative success, a rigorous, formal understand- 
ing of the level of fluctuations included by the "classical field 
method," and its relation to more traditional CF methods, is 
still lacking. 

In this paper, we establish that, in the case of an interacting 
Bose gas in a one-dimensional ring geometry, the "classical 
field method" approaches the same equilibrium as the more 
traditional stochastic and molecular dynamics CF methods, 
and that this equilibrium agrees with an exact evaluation of 
the partition function in the high-temperature limit. First, we 
discuss the stochastic dynamics approach which, in the flavor 
of the time-dependent Ginzburg-Landau equation, relaxes the 
field configuration to the minimum of the free energy while 
accounting for thermal fluctuations about this minimum. Sec- 
ond, we evaluate the microcanonical dynamics method which 
resembles a molecular dynamics approach based on Hamil- 



ton's equations of motion. Third, we look at the "classical 
field method" which is based on the non-linear Schrodinger 
equation (NLSE). 

In the case of the one-dimensional Bose gas on a ring, that 
we consider here, there is a simple way to evalutate the parti- 
tion function at finite temperatures known as the transfer inte- 
gral method. This solution serves us as a reference to which 
we can compare the properties of the steady states which the 
simulations reach. We find that the thermodynamic proper- 
ties, such as distribution and correlation functions, of all three 
simulations, agree very well with the exact solution given by 
the transfer integral. 



II. THE TRANSFER INTEGRAL SOLUTION 

We begin by evaluating the high-temperature limit of the 
partition function for an interacting Bose gas on a one- 
dimensional ring, exactly, by means of the transfer integral 
method [6]. This will serve as a reference by which we may 
judge the various CF methods discussed in this paper. 

All information about the equilibrium properties of a many- 
body system may be extracted from its partition function Z = 
J D(<fi, <fi*)e~~ s , where the action S[</), cf>*] reads 



S[<t>, 4>*] = J drjdx (<f>*d T + H[<j>, ~ ^[0, 0*]) , 

(1) 

fi is the chemical potential, and (3=1 /T denotes the recipro- 
cal temperature (we set ks = 1). 

In the high temperature limit, j3 — ► 0, the integration over 
imaginary time becomes quenched, and paths not on the t = 
axis can be neglected, thus one arrives at what we will refer to 
as the "classical partition function" 



(2) 



where we have introduced F = J dx (H — fiN) the free en- 
ergy. Although simpler than the full quantum-mechanical par- 
tition function, the classical partition function of Eq. © is, 
in general, a very difficult quantity to calculate. This is es- 
pecially true near a second-order phase transition where all 
orders of classical fluctuations may contribute. 
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A Bose gas with contact interactions can be described by 

F= fdx [0* (-fr 2 V 2 /(2m) - fjL)(j> + «?|0| 4 /2] , (3) 



where g is the interaction constant and m is the atomic mass. 
Its free energy is of the general form 



F = J dx {a\(j>\ 



2 + bU\ A + c\Vcb\ 2 ) 



(4) 



Defining a length scale £o = y c / a = h/y / 2~mjl which in- 
duces an energy scale e = ^ 2 /(2™£ 2 ) we can rewrite the 
free energy in dimensionless quantities 

PF = pJdx (a|^| 2 + 6|^| 4 + |V^| 2 ) . (5) 

In the following we will leave off the tilde and all symbols 
refer to dimensionless quantities. 

In one spatial dimension, the transfer integral method re- 
lates the classical partition function of Eq. (O, for a free en- 
ergy of the form of Eq. @, to the Schrodinger equation of a 
single particle in an anharmonic oscillator potential 



4/3 2 



V 2 U + a\u\ z + 6|u| 4 tp n (u) = E n ip n {u). (6) 



In the appendix we show that thermodynamic properties, such 
asthedistributionfunctionoffieldamplitudesP(|0| = it) and 
the two-point correlation function (</)(•)* <fi(»+r)), can readily 
be calculated from the eigenfunctions ip n (u) and eigenvalues 
E n of the single-particle Schrodigner equation, Eq. ©, by 
means of 



and 



Pi 



+ r)) 



u) = 2u\ip (u)\ :: 



E 



du <fin( u ) U <Po( u ) 



-P{\r\/to){E n -Eo) 



(7) 



(8) 



The Schrodinger equation of Eq. © can easily be solved nu- 
merically by exact diagonalization. This approach provides 
a powerful way to evaluate the high-temperature limit of the 
partition function for an interacting system exactly which we 
can compare to our numerical simulations. 



III. STOCHASTIC DYNAMICS 

An alternative method to evaluate the partition function of 
Eq. (0 is to allow the fields to evolve in such a way that they 
relax to the correct equilibrium distribution and to dynami- 
cally sample the relevant observables. This method of deter- 
mining the equilibrium quantities was originally proposed by 
Parisi and Wu |7|], while a related analysis has been used to 
model the behavior of a Fermi gas close to the BCS transition 
point d. 



The fictional "dynamics" are determined by a diffusive term 
which relaxes the field configuration to the minimum of the 
free energy surface while a white noise term randomly drives 
these fields. The result is a Langevin equation for the dynam- 
ics of the fields 

d t 4> = -{acj> + 2b\(j>\ 2 (t>-V 2 <f>) +£, (9) 

where the noise correlations are those of white noise 

(e(x,t)Z(x',t')) = ^8(x-x')5(t-t'). (10) 

Note that this is a grand-canonical method, so the inverse tem- 
perature (3 appears explicitly. 



A. The non-interacting case 

Let us first look at the non-interacting case (b — 0) in which 
Eq. (O takes on the form of an Ohnstein-Uhlenbeck process 
and is analytically solvable |9J]. 

For long times, t — > oo, the distribution function of field 
amplitudes P(Re(<fi) — u) approaches 



P(Re(» =u) = 




-f3\/2au 2 



and the spatial correlation function becomes 



1 



(<X.)>(. + r))= . , 
which in momentum space reads 

1 



-v^(|r|/{ ) 



mok) 2 + 2a) : 



(ID 



(12) 



(13) 



expressing the equipartition of energy between the classical 
modes with the free dispersion relation E(k) = k 2 . 

Using a second-order Runge-Kutta method, we propagate 
the initial state <f>(x, 0) = with the Langevin equation ©, 
where the stochastic noise ( TTOb is interpreted in the Ito sense 
|9fl. After the observables, i.e. the distribution and correlation 
functions, have settled to a steady state, we calculate their val- 
ues by averaging over a set of trajectories and compare them to 
the analytic expressions in Eqs. ( fTTI ). ( fT2l and ( fT3l ). We used 
a uniform spatial grid with step size Ax = 0.2, an interval of 
length L = 20, a time step At = 0.05Ax and the parameters 
P = 2, a = 1/2, b = 0, similar to the work in Oil]. We av- 
eraged over 1000 trajectories at a final time of tf = 20, when 
the system has reached a steady state. 

In Fig. Q] we show the results of our numerical simula- 
tions. The distribution function for the real part of the field 
P(Re(</>) = u) and the distribution function for the modulus 
of the complex field amplitude P(\<f>\ = it), as well as the cor- 
relation functions in real space ((f)(9)* (/)(• + r)) andits Fourier 
transform (\<t>k\ 2 ) are presented aside with the analytical solu- 
tions (Qj), CE3 and ( fT3b . 

The equilibrium properties of the numerical simulation of 
Eq. (O agree extremely well with the analytic solutions. We 
also used this baseline test to check the transfer integral ap- 
proach, which overlaps exactly with the analytical solutions. 
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FIG. 1: (Color online) Langevin dynamics in the non-interacting 
case (/3 — 2, a — 1/2, b = 0). (a) Radial distribution function 
P(\4>\ ~ u) and (inset in (a)) Distribution function P(Re(<j!>) = it), 
(b) Correlation function in real space {</)(•)* 4>(» + r)) and (inset 
in (b)) in momentum space (\<f>k\ 2 )- Numerical results (points) from 
the steady state of the Langevin equation ^ and analytic results from 
{ED, G2J and {13} (solid line). 
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FIG. 2: (Color online) Langevin dynamics in the interacting case 
(J3 = 2, a = -1/2, 6 = 1/4). (a) Radial distribution function 
P(\4>\ = u ) an d (inset in (a)) Distribution function P(Re(0) = u), 
(b) Correlation function in real space (<^(»)*<^(» + r)) and (inset in 
(b)) in momentum space (|</>fc| 2 ). Numerical results (points) from 
the steady state of the Langevin equation ^ and exact solution from 
the transfer integral method (|6l( (solid line). 



B. The interacting case 

With the transfer integral method we can extend our stud- 
ies to the non-linear case (6 ^ 0), where results cannot be 
obtained analytically. In Fig. [2] we plot the distribution and 
correlation functions for (3 = 2, a = — 1/2 and b = 1/4. All 
other parameters were the same as before. We see that the 
stochastic simulation correctly reproduces the exact thermo- 
dynamic results for this interacting system in the limit of long 
times where the system has settled to a steady state. 

Fixing all parameters, we explore the equilibration process 
as we decrease the temperature, i.e. increase (3. In Fig. |3]dis- 
tribution and correlation functions are shown for (3 = 6. We 
find that the time to reach the steady state increases by over an 
order of magnitude, from tj = 20 to tf = 250. In the real part 
of the distribution function we can identify the minima of the 
Mexican hat potential V(x) = ax 2 +bx 4 to which the dynam- 
ics is limited for small temperatures. The distribution function 



P(\4>\ = u) still agrees well with the exact solution (0. While 
the simulation equilibrates to the correct density (!</>(•) | 2 }, its 
spatial correlation function (</>(•)*</!>(• + r)) deviates for large 
spatial scales from the transfer integral as given in Eq. ([8]). 

To explain this behavior we recall that the long-range prop- 
erties of the correlation function are determined by the two 
smallest eigenvalues of the transfer integral: Eq < E\. We 
therefore define a correlation length l c by 



Co /3(£i-£ )' 

In Fig.|4]the correlation length l c is calculated from the trans- 
fer integral of Eq. (|6]l as a function of the inverse temperature 
(3. The correlation length l c becomes comparable to the sys- 
tem size 3l c w L/2 for f3 > 4. At this temperature coherence 
is established across the whole system which is not accounted 
for by our expression for the correlation function as given by 
Eq. I©. Instead, we must properly account for the finite size 
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FIG. 3: (Color online) Langevin dynamics in the interacting case 
(J3 = 6, a = -1/2, 6 = 1/4). (a) Radial distribution function 
P{\(f>\ = u) and (inset in (a)) Distribution function P(Re(<j>) = it). 
Numerical results (points) from the steady state of the Langevin 
equation l|9) and exact solution from the transfer integral method ((6) 
(solid line), (b) Correlation function in real space (<£(•)*<£(• + r)) 
aside with numerical results (points). The transfer integral solution 
is presented including finite-size corrections dA. lOfc (solid line) and 
in the thermodynamic limit l[8j (dashed line). 



of the ring in our derivation of the correlation function (which 
we show in the appendix). Using Eq. (IA. 10b . we plot in Fig. [3] 
the exact result for the correlation function of this finite sys- 
tem. The simulation data agrees very well with this version 
of the correlation function. Moreover, we note that the long 
equilibration times observed in the simulations at lower tem- 
peratures are explained by the fact that it takes many collisions 
to establish coherence between the particles, the typical 'slow- 
ing down' phenomenon close to a phase transition. 



IV. MOLECULAR DYNAMICS 

One of the disadvantages of the stochastic method we have 
just presented is that it is difficult to generalize to an inhomo- 
geneous situation. Here we present an alternative CF method 
which (1) may easily accommodate an inhomogeneous poten- 
tial, (2) is a microcanonical method so may be more directly 




FIG. 4: (Color online) Correlation length l c as a function of the in- 
verse temperature (5 for a complex field (solid line). Additionally, we 
plot the correlation length l c for a real field (dashed line). For real 
fields, the correlation length increases much more dramatically with 
decreasing temperature and, therefore, equilibration at temperatures 
/3 > 6 requires very long interaction times. 



compared to the classical field method, and for which (3) an 
effective temperature of the interacting system may be unam- 
biguously determined because the momentum field is decou- 
pled from the density field. 

Starting from the partition function of Eq. (01, one calcu- 
lates observables by the relation 



(0> = 



fD(cj>,ct>*)e-S 



(15) 



We introduce the canonical momentum field n = d t (f>, which 
is independent of the density field (j>, by multiplying the parti- 
tion function of Eq. (0 by unity 



D(ir,Tr*)e- l7Tl2 oc 1 



and arrive at the partition function 



/? / dxH[ir, 



(16) 



Z oc / D(Tr,n*) I L>O,0*)exp 

J (17) 

where H [ir, <j>] may be treated as a Hamiltonian with the action 
S[4>] playing the role of a potential 

Hfcfl = \tt\ 2 + S^}. (18) 
We may readily derive Hamilton's equations of motion 

7T = -(a4> + 2b\(j)\ 2 (t>- V 2 0) , (19) 



IT, 



and combine them to obtain 



d?(j> = V 2 cf>- (Hp- 2b\(f)\ 2 (j>. 



(20) 



(21) 
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FIG. 5: (Color online) Molecular dynamics. Kinetic energy per par- 
ticle (|7r| 2 ) versus time t and (inset) correlation function in Fourier 
space in the steady state (\(j>k\ 2 }- Simulation results (points) and ex- 
act solution (solid line). 



Using a standard leap-frog integration routine we solve 
Eq. (BTI ) numerically for the non-equilibrium initial condition 

iV EX 

0(x,O) = 0, 7T(>,0) = y £^Ae 2 - Ki{ - kx l L -^\ (22) 

k=l 

where k is an integer and are random numbers drawn from 
a uniform distribution of < < 2tt. The initial energy 
is then simply E = N CX LA 2 Hill . For our simulation, we 
choose, similar to the last section, iV ex = 4, E/N = 2, a = 
-1/2, and b = 1/4 as well as L = 20, Ax = 0.2, and At = 
0.05 Ax. 

In Fig. [5] we show the average kinetic energy per degree of 
freedom (number of grid points in the simulation) as a func- 
tion of time. The system comes to a steady state after mixing 
for 2500 time units. In the steady state, we compute the cor- 
relation function (\<t>k\ 2 ) by replacing ensemble averages with 
time averages over 500 time units and, additionally, over 10 
trajectories, to further reduce fluctuations. We may now de- 
termine an effective temperature T e ff of the interacting system 
by computing the time-averages of the kinetic energy (|7r| 2 ) 

T eff = (|tt| 2 ) Ax/2 = 0.47. (23) 

In the inset of Fig. [5] we plot the correlation function from 
our simulations aside with the exact solution obtained from 
the transfer integral approach for (3—1 /T e g. The agreement 
between the properties of the steady state and the exact solu- 
tion, as well as with the simulations of the previous section, is 
very good. 

We would like to stress that since the tt field is decoupled 
from the density field <j>, a temperature may be reliably deter- 
mined from these microcanonical simulations. This, unfortu- 
nately, is not the case in the following microcanonical method. 



FIG. 6: Schematic diagram showing the classical and incoherent re- 
gions of the single particle spectrum for a harmonically trapped Bose 
gas. 

V. "CLASSICAL FIELD METHOD" 

Let us now approach the interacting Bose gas from a dif- 
ferent perspective. To begin, we will start from the real-time 
quantum action for an interacting Bose gas in the microcanon- 
ical ensemble 

S[4>, <t>*] = -<j>*idt<t> + b\<f>\ 4 + |V</>| 2 . (24) 

Minimizing the action, which is equivalent to making a 
saddle-point approximation to the full quantum mechanical 
kernel, we arrive at 

O=^ = -z9 t + 26|0| 2 0-V 2 0, (25) 

i.e. the usual non-linear Schrodinger equation (NLSE). At this 
level of approximation the dynamics of the fields are, roughly 
speaking, restricted to classical paths, and one might think 
that a microcanonical ensemble of particles, governed by such 
an equation of motion, displays equilibrium properties deter- 
mined by the classical partition function of Eq. ©. 

Note, the NLSE equation is usually thought to account only 
for the lowest energy mode of a Bose gas, but this restriction 
is nowhere imposed in the derivation above. The validity of 
the saddle-point approximation, however, will depend upon 
the extent to which each mode can be treated classically. In 
Fig. [6] we depict the energy levels of a harmonically trapped 
Bose gas. In a finite region, from the ground state up to some 
cutoff energy, all the modes are multiply occupied. In this 
region, quantum fluctuations are negligible and we may treat 
these particles as a set of coherent or classical fields. For more 
details as to the subtleties of this method, such as how to in- 
corporate the incoherent portion of the spectrum, we refer the 
reader to the literature [12, 

El [II El. 

We start our simulation from the non-equilibrium state 

2-ninx 

4>{x)= -7r e ^ (26) 
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FIG. 7: (Color online) Classical Field Method. Steady-state momen- 
tum distribution at C = 2b — 1000 (points) and the new steady- 
state after the change in the interaction constant to C — 2b = 5000 
(circles). Both are fitted to the exact solution given by the transfer 
integral {6]l (solid and dashed line). 



where 2N ex + 1 is the number of excited modes, L the sys- 
tem length and £ n are random numbers drawn from a uniform 
distribution in the interval [0, 2ir}. We propagate this initial 
condition with the NLSE (l25T > using a standard semi-spectral 
fourth-order Runge-Kutta routine with L = 20, Ax = 0.2 
and At = 10~ 4 and choose the non-linear coefficient to be 
C = 2b = 1000. After the NLSE d25]l reaches a steady state 
lfl2Tl . we take time averages over a period of 1000 time units 
to calculate the momentum distribution (|0j,| 2 ). To compare 
with the exact transfer integral solution, we first adjust the 
chemical potential a to match the particle number. Unfortu- 
nately, we do not know of a simple way to directly extract the 
temperature from these simulations when the interactions are 
large, as in the case we are currently concerned with. How- 
ever, we may avoid this issue by fitting the curves to the exact 
solution, treating the the temperature j3 as a free parameter. 

In Fig. |7] we show the momentum distribution of the steady 
state, to which the classical field method equilibrates at strong 
non-linearities, and compare it to an exact solution for j3 — 
5.29, a = —63.1 and b = 1000/2. The agreement is excel- 
lent which strongly supports the conclusion that the classical 
field method is able to correctly sample the classical partition 
function. 

If, after equilibrium has been reached, we suddenly change 
the interaction constant to a different value, C — 2b = 5000, 
we find that the momentum distribution changes to the new 
equilibrium value which can be fitted with the exact solution 
j3 = 2.55, a = -313.3 and b = 5000/2. We note that the fit 
does not reproduce the tails of the distribution exactly, how- 
ever, the overall fit is still excellent. 

We believe that these observations are a strong motivation 
for using the "classical field method", since it is based upon 
the NLSE, to describe dynamics of cold atomic gases at finite 
temperatures. This subject, however, is beyond the scope of 
the current paper. 



VI. CONCLUSION AND OUTLOOK 

We have compared three different CF approaches to dynam- 
ically sample the classical partition function to an exact trans- 
fer integral solution. First, a stochastic dynamical method, 
based upon the evolution of a Langevin equation, was shown 
to contain the correct thermal fluctuations needed to recover 
the exact distribution and correlation functions of the full clas- 
sical theory. Second, we presented the microcanonical method 
of molecular dynamics, which was again able to reproduce 
the exact one-dimensional solution. Both of these dynami- 
cal schemes have either a temperature built-in or one which 
may easily be determined - although the latter is more eas- 
ily adapted to an inhomogeneous system. Third and final, we 
have evaluated an alternative CF method, the "classical field 
method", and shown that its steady state is given by the ex- 
act result for the partition function of a Bose gas in the high- 
temperature limit. 

Furthermore, when we changed the interaction strength of 
our "classical field method" simulations, the system was seen 
to dynamically evolve to another equilibrium point even in 
this low-dimensional system. This observation provides a 
strong motivation that the classical field method may be used 
to study finite-temperature dynamics. It would be especially 
interesting to apply this formalism to experiments in ID ring 
geometries with Josephson junctions which can be realized 
e.g. by perturbing the homogeneous ring by a potential bar- 
rier. In this case, the molecular dynamics approach should be 
valuable as it can easily be extended to this inhomogeneous 
situation and provides an effective temperature for the inter- 
acting system, while the classical field method offers a power- 
ful means to study the dynamical relaxation process. Another 
line of research might be the study of heat conduction in trans- 
port experiments on strongly-interacting systems. 



APPENDIX: DERIVATION OF THE TRANSFER INTEGRAL 
SOLUTION 

In this appendix we review the evaluation of the classical 
partition function, Eq. {2), by means of the transfer integral 
Og]. Assuming a free energy of the form of Eq. (0, we next 
discretize the one-dimensional domain [0, L] with N equally 
spaced grid points and apply periodic boundary conditions 
<A/v+i = <pi> turning the domain into a ring of length L. The 
partition function then takes the form 



N 



-|8(Aa:/£o)/(&+i,<M 



where 



f(<t>i+u<k) = a \4>i+i\ 2 + fr|<A;+i| 4 + 



Ax 



By introducing the following generalized notation 



,-/3/(^-,^) 



ri 



(A.l) 



(A.2) 



(A.3) 
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the partition function can be rewritten as 



Z = tT 



T 



N 



(A.4) 



If we can find the eigenfunctions \ip n ) and eigenvalues E n of 
the transfer operator T, the partition function can be exactly 
evaluated through the relation 



n=0 



-P(L/H )E n 



(A.5) 



The eigenvalue problem T \ <p n ) = \ n \tpn) reads explicitly 



e-Mti+uM^fa) = e -^^„(^ +1 ), (A.6) 



and can be restated by expanding ip n (4>i) in a Taylor series 
about tp n (4> i+1 ) as 



e-< iB <p n = e-f» s " Vn , 



(A.7) 



where H is the Hamiltonian of a single particle in an anhar- 
monic oscillator potential 



H = 



In this way we have reduced the problem of calculating the 
classical partition function (O to finding the eigenvalues and 
eigenfunctions of a simple, one-particle, time-independent 
Schrodinger equation H(f n = E n (p n . 

The spatial correlation function (</>(•)*(/>(• + 7-)) can be de- 
termined as follows 



(<f>(.)*<t>(.+r)) = Z~ l J DfaPMoytWe 
= Z- 1 tr U(0)f r (j){r)f N - 



-/3F 

(A.9) 



Using the spectral representation of T = \(p n )^n( t fin\< we ob- 
tain 



, JV-r 



(0(.)>(. + r)) 

E l3 \(^\^)\ 2 (e^Y(e~^) N - r 



(<Pi\0(r)\<fj) 



(A. 10) 



where we have used the spatial homogeneity of the system to 
arrive at the last line. 

In the thermodynamic limit, N — > oo, and when the system 
size is large compared to the coherence length l c <C L, we can 
simplify this formula to 



+ r)) = |(^|0l^o>| 2 e ^ (fi ^ £o)(r/ « o) , 

i=l 

(A. 11) 



(A. 8) where the matrix element has to be evaluated through 



((pi\<t)\(po) = du<p*(u)ucp (u). 
Jq 
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